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ABSTRACT 

We investigate the physical mechanisms of tidal heating and satellite disruption in 
cold dark matter host haloes using N-body simulations based on cosmological initial 
conditions. We show the importance of resonant shocks and resonant torques with the 
host halo to satellite heating. A resonant shock (torque) couples the radial (tangential) 
motion of a satellite in its orbit to its phase space. For a satellite on a circular orbit, 
an ILR-like resonance dominates the heating and this heating results in continuous 
satellite mass loss. We estimate the requirements for simulations to achieve these 
dynamics using perturbation theory. Both resonant shocks and resonant torques affect 
satellites on eccentric orbits. We demonstrate that satellite mass loss is an outside-in 
process in energy space; a satellite's stars and gas are thus protected by their own 
halo against tidal stripping. We simulate the evolution of a halo similar to the Large 
Magellanic Cloud (LMC) in our Galactic dark matter halo and conclude that the LMC 
stars have not yet been stripped. Finally, we present a simple algorithm for estimating 
the evolution of satellite mass that includes both shock heating and resonant torques. 

Key words: galaxies: evolution — galaxies: interaction — galaxies: haloes — galax- 
ies: kinematics and dynamics — method: numerical — method: N-body simulation 



1 INTRODUCTION 

Physical processes affecting satellite galaxy evolution in 
their host haloes are an important component of galaxy for- 
mation in the cold dark matter (CDM) cosmogony as galax- 
ies are built up from the assembly of small structures. This 
assembly includes the process of satellite galaxies merging 
with their host galaxies. Moreover, recent CDM cosmolog- 
ical simula tions predict the ex i stence of a large n umber of 
subhaloes (|Ghigna et al.lll998l ; Iklypin et all 119991 ). Conse- 
quently, understanding the detailed physical processes af- 
fecting satellite evolution are key ingredients to understand- 
ing galaxy formation in the CDM cosmogony. 

Several basic questions about satellite halo evolution 
remain. First, how is the satellite stripped? In other words, 
what parts of the initial mass distribution might persist to 
the present day? Second, what is the rate of satellite halo 
disruption? Third, how does a satellite's internal structure 
evolve? In a satellite-galaxy merger, the stars and gas of 
the satellite galaxy can be stripped and become halo stars 
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and halo gas (e.g. iQuilis et all l2000l : iBullock et all 1200 ll ). 

The interaction between the satellite galaxy and its host 
during the course o f the merger resul ts in evolution of the 
satellite galaxy fe.g. lMoore et al1ll996l ). The remaining com- 
ponents of the satellite galaxy merge with the host galaxy 
and this merging causes the host galaxy to gain mass (e.g. 



and this merging c 
iMurali et all 12002) 



Current cosmological simulations can only provide sta- 
tistical properties of subhaloes since even in the high- 



est resolution cosmological simulations llGhigna et al 
| De Lucia et alj|2004 IDiemand et"al]|2004l: IGao et al 



2000; 



2004 



l Oguri fc Led 120041 : IDiemand et all 120071 : ISpringel et ail 
I200SI ) the detailed physical processes of individual sub- 



halo evolution have not been accurately studied, owing 
to limited resolution. To investigate the physical pro- 
cesses in detail, we perform high resolution idealised 
simulations with cosmologically motivated initial condi- 
tions instead of using cosmological simulations. In these 
idealised simulations, a live satellite orbits in a static 
host halo. Although too simplified to reproduce a satel- 
lite's evolution in realistic detail, several authors have 
used similar non-cosmological simulations to study satel- 
lite disruption with alternative simula t ion methods (e.g. 



Havashi et all 120051: iKazantzidis et all 120041: iRead et al.1 



2006: 



iBovlan-Kolchin fc Mal l2007l ; IPenarrubia et all 12008) 
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These studies have demonstrated important features of the 
satellite evolution such as mass loss history and density pro- 
file evolution. In this paper we use the higher resolution 
simulation to investigate satellite disruption particularly fo- 
cusing on the detailed physical processes affecting satellite 
evolution such as resonant dynamics. 

When a time-dependent force acts on a bound system 
such as a galaxy or a dark matter halo, resonant interac- 
tions play an important role in the system evolution. Recent 
studies of resonant dynamics in galaxy evolution claim that 
high resolution simulations are required to accurately repro- 
duce these resonant effect s (|Weinberg fc Ka tz 2007a. bj) and 
IWeinberg fc Katzl (|2007al ^ provide a procedure to determine 
minimum particle number guidelines. Since our idealised, 
high resolution simulations are designed to satisfy these par- 
ticle number guidelines, they allow us to investigate the role 
of resonant dynamics in satellite disruptions. 

Resonant interactions couple a time-dependent perturb- 
ing force with orbits in the system. The frequency spec- 
trum of the time-dependent perturbing force characterises 
the interaction. For a satellite orbiting in its host halo, the 
satellite's orbital frequencies, and possibly the rate of or- 
bital decay, determine the time-dependence of the external 
force. A general, eccentric satellite orbit in a spherical halo 
has both a radial and azimuthal frequency, making the res- 
onant coupling for an eccentric orbit complex. Empirically, 
we may characterise the overall effects of the interaction as 
a resonant shock and a resonant torque. A resonant shock 
represents coupling with the radial orbital frequency and a 
resonant torque represents coupling with the azimuthal or- 
bital frequency. 

A resonant shock is a generalisation of the standard im- 
pulsive, gravitational shock. During a resonant shock, some 
orbits within the satellite gain energy through resonant cou- 
pling even though they are not in the impulsive limit, i.e. the 
time scale of the perturbation near pericentre is much longer 
than the internal orbital time scale. A resonant torque cou- 
ples the rotation in the external potential to orbits in the 
satellite. For a resonant torque, the magnitude of the exter- 
nal potential does not have to change; a change in the posi- 
tion angle of the satellite frame relative to the centre of the 
halo is sufficient to produce a torque. Resonant shocks have 
been previously considered and included in the impulsive 
approximation a s an adiabatic correction (e. g. ISpitzer|[l987l ; 
Weinberg 1994allbi: iGnedin fc Qstr ikcr 1999). However, res- 
onant torques have not been similarly considered in satellite 
evolution studies although they have been extensively inves- 
tigated in the dynamics of barred galaxies. In this study, we 
will carefully investigate these resonant effects on satellite 
evolution using this distinction. 

Globular c luster evolution in a host galaxy is well es - 
tablished (e.g. ISpitzed 1 19871 ; IChernoff fc Weinberg! Il990t ). 
A globular cluster experiences both tidal truncation and 
heating by both compressive gravitational shocks and tidal 
shocks. Because satellite halo evolution in a host halo is 
similar to globular cluster evolution in a host galaxy, many 
galaxy formation studies employ simple analytic formulae 
taken from these globular cluster evolution studies to es- 
timate satellite galaxy evolution. However, unlike globu- 
lar clusters, the satellite-host mass ratio is not vanishingly 
small. This breaks the spatial symmetry in mass loss, as 



shown in IChoi et al.ll|2007h . and changes the relative impor- 
tance of resonant coupling. 

In this paper, we present numerical simulation results of 
satellite galaxy disruption. In SJ2] we present an overview of 
our numerical techniques: the N-body simulation code, the 
generation of initial conditions, and the relevant perturba- 
tion theory. In §3] we present the results of a circular orbit 
simulation. We show that resonant torque effects result in 
significant satellite mass loss. In §31 we present the results of 
eccentric orbit simulations. We show that satellite heating by 
gravitational shocks at pericentre, which also includes inter- 
nal structure evolution, is the dominant process responsible 
for disrupting the satellite. In fj5l we show that the process of 
satellite stripping is an outside-in process in satellite energy 
space. Using this finding, we suggest an explanation for the 
'missing' LMC stellar stream. We also discuss the evolution 
of the satellite density profile. In §6] we provide an improved 
analytic estimate for satellite mass loss, and we summarise 

in m 



2 NUMERICAL METHODS 

Our N-body simulations use a three-dimensional self- 
consistent field code (SCF, also known as an expan- 
sion c ode; |ciutton-Brocklll972l , Il973l ; iHernquist fc Ostrikerl 
1 19921 ; IWeinberg||l999h . The expansion code calculates bi- 
orthogonal basis sets of density-potential pairs and com- 
putes the gravitational potential of the system using these 
basis sets. There are two reasons that an expansion code 
is an appropriate potential solver for our study. First, the 
expansion basis can be chosen to follow the structure over 
an interesting range of scales and simultaneously suppress 
small-scale noise. In contrast, the noise from two-body scat- 
tering can arise at all scales in direct-summation, tree, and 
mesh based codes. Small-scale noise can also give rise to 
diffusion in conserved quantities, which can then lead to 
unphysical outco mes particularly for stud ies of long-term 
galaxy evolution l|Weinberg fc Kat z 2007a b|). Second, the 
expansion code is computationally efficient; the computa- 
tional time increases only linearly with particle number and 
with modest overhead. Hence, the expansion code permits 
the use of a much larger number of particles than most other 
codes for the same computational cost. 

The expansion code is not adaptive. The largest effi- 
ciency obtains when the basis resembles the galaxy. Oth- 
erwise, the expansion code requires a large number of basis 
function pairs, which introduces small-scale noise and results 
in a greater computational overhead. These constr aints have 
been minimised by so me recent improvements l|Weinberel 
1999]; [Ch oi et al.l I2007T ). First, employing a numerical so- 
lution of the Sturm-Liouville equation, the initial galaxy 
model can be used as the zeroth-order basis function for 
the expansion code. Then, the expansion code requires only 
a modest number terms to accurately compute the poten- 
tial. Second, we accurately trace the density-potential cen- 
tre of the expan sion during the course of a simulation (see 
IChoi et alj |2007| . for details). Third, to reduce truncation 
error, we separately track the centre of mass motion of the 
satellite and the relative motion of the satellite centre. With 
these improvements, the expansion simulations are now able 
to achieve a sufficiently high central resolution with suffi- 
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ciently low noise to allow the investigation of the detailed 
dynamical processes important for satellite disruption. 

All the halo models used in our simulations are 
based on the cosmol ogically-motivated universal CDM halo 
(|Navarro et al.lll997l , hereafter NFW), p(r) oc r' 1 (r + r s )~ 2 , 
where r s is a scale length characterised by the concentration 
parameter c = R V i r /r s and R V i r is the virial radius of the 
halo. We represent the host halo potential by a concentra- 
tion c = 15, static NFW halo. The NFW profile has infinite 
extent but a real satellite is tidally truncated. Furthermore, 
the host halo tidal field will affect a satellite halo even before 
the satellite halo passes within the host halo's virial radius, 
but it is computationally expensive to simulate the evolution 
of a satellite from such large radii. Since the objective of our 
study is only to understand the physical processes respon- 
sible for satellite disruption, it is not necessary to simulate 
the satellite's entire evolutionary history. Therefore, we be- 
gin the simulation with the satellite in the host halo on the 
desired orbit and include the host halo tidal field when we 
generate the satellite's initial conditions. 

Our simulations ignore the effects of dynamical friction 
and the subsequent reaction of the host halo since the satel- 
lite masses of interest are oft en much smaller than that o f 
the host halo. According to (|Bovlan-Kolchin fc Mai l2007h . 
when the satellite-to-host mass ratio is smaller than 0.1, the 
dynamical friction hardly makes the satellite decay within 
half a Hubble time. The gravitational back-reaction of the 
host halo has a complex structure with an amplitude pro- 
portional to mass ratio. In many cases, the halo reaction is 
dominated by barycentric motion of the halo center induced 
by the satellite perturbation. We will quantitatively discuss 
these dynamics in a later paper. Since our main objective in 
this paper is an understanding of the physical processes of 
satellite disruption, we do not expect omission of the halo 
response to affect our conclusions. 

The details of the procedure that we use to generate a 
satellite's initial conditions are as follows. We start with an 
initial NFW satellite halo model. Since the maximum cir- 
cular velocity is a better measure of a satellite's initial size 
than its mass, owing to ongoing mass loss, we characterise 
the satellite's size by its circular velocity V = \J M V i r / R v ir ■ 
We choose the initial satellite halo model such that the circu- 
lar velocity of the satellite V sa t = \Vhost, which means that 
the initial virial mass of the satellite is 0.0046. Throughout 
this paper we use system units unless otherwise specified 
with G = 1, M virthost = 1, and R vir , host = 1. If scaled to 
the Milky Way, this satellite halo corresp onds to the size 
of th e Sagittarius dwarf dark matter halo (Maiewsk i et al.l 
12004! ) . We then determine the truncation radius, for which 
we discuss several choices below, assuming that the satel- 
lite is on a circular orbit with R = 0ARh os t,vir- We call 
this fixed distance the tidal distance. As shown in Fig. [TJi, 
the truncated model is identical to the initial NFW model 
but is chopped at the truncation radius. After truncating 
the initial satellite model at this radius, we perform an Ed- 
dington i nversion to compute the co rresponding distribution 
function (|Binnev fc Trem ainc 1987, Chapter 4), as shown in 
Fig- [It- We calculate the final satellite density profile by in- 
tegrating this distribution function over velocity. In contrast 
to the truncated model, the final satellite model does not 
show a sudden drop at the outer edge of its density profile 
(see Fig. [TJ;). Owing to this smooth outer profile, the final 
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Figure 1. The effect of our truncation procedure on a satellite's 
initial NFW profile, (a) The density profiles of the initial NFW 
model and the truncated model, (b) The distribution function 
computed using Eddington inversion, (c) The density profiles of 
the truncated model and the final satellite model, (d) The circular 
velocity profiles of the initial NFW model and the final satellite 
model. The truncation radius in this figure is x e and the satellite 
is on a circular orbit at R = 0AR v i r . 



satellite model is closer to equilibrium than the truncated 
model. We use this final satellite model as our initial satel- 
lite halo model. Fig. QJl shows the circular velocity profiles 
of the initial and final satellite model; even though the final 
satellite model is considerably truncated compared to the 
initial model and has lost 2/3 of its initial mass, the maxi- 
mum circular velocities of both models are nearly the same. 
We realise the satellite phase space using an acceptance- 
rejection method. All realisations are assumed to have an 
isotropic velocity distribution and are made of 10 6 equal- 
mass particles. 

The Roche limit, th e critical-poin t radius x e , and the 
tidal cut-off radius (rt) |Spitzerl [l987. Chapter 5) are all 
natural choices for the truncation radius. The Roche limit 
is the radius where the satellite density and the host halo 
density are the same. The critical-point radius, x e , defines 
the unstable zero-velocity equilibrium; that is, the unstable 
equilibrium in the effective potential d$ e // = 0. For ease 
in computing initial conditions and idealised tests in later 
sections, we define a sphericised centrifugal potential term 
of the form 

1 2 
-il ar 
2 

with a £ [0, 1] (see Appendix |A1 for details). For initial con- 
ditions, we set a = 1 to make the satellite smaller and hence 
to minimise satellite mass loss owing to initial adjustments 
when it is first placed in the external potential. The tidal 
cut-of£| is the radius where rt = 2x e /3. This is the radius of 

We follow Spitzer (1987) in defining the tidal cut-off radius 



© 0000 RAS, MNRAS 000, 000-000 



4 Choi, Weinberg, & Katz 





Figure 2. Density profiles for three satellite initial conditions 
truncated at three different radii: the Roche limit (red), x e 
(green), and rt (blue). 



the critical potential surface perpendicular to the direction 
toward the host halo centre so it is the maximum radius 
that a spherical satellite halo can have and not extend be- 
yond the critical potential surface in any direction. We plot 
the density profiles of the truncated satellite halo models in 
Fig. [2] To obtain the final density profiles for the satellite 
halo models truncated at x e and r t we apply the Eddington 
inversion procedure. However, we do not apply this proce- 
dure when we truncate at the Roche limit since it would 
change the outer density, making it differ from that of the 
host halo, and it would no longer be the Roche limit. Fig. 
[2] shows that at the given tidal distance of 0.4, the satellite 
truncated at the Roche limit is the largest and the satellite 
truncated at rt is the smallest. 

Fig. [3] shows the evolution of mass from circular orbit 
simulations for the three satellites shown in Fig. [2] (see ij3] 
for a detailed discussion). All three satellites show continu- 
ous mass loss. The amount of mass loss correlates with the 
satellite size: the satellite truncated at the Roche limit loses 
the largest mass fraction and the satellite truncated at r t 
loses the least mass. We choose x e for our fiducial trunca- 
tion radius on physical grounds. x e is the transition point 
beyond which satellite particles on a circular orbit become 
unbound. At the Roche limit radius for an NFW halo, the 
host halo potential dominates the gravitational potential of 
the satellite. In contrast, using the tidal cut-off radius, r t , 
is a sensible choice for globular clusters, which have already 
orbited many times around a galaxy and have already been 
severely truncated, but using r t could underestimate a satel- 
lite halo's size since it has made only a few complete orbits. 
Hence, we restrict our study to initial satellite haloes that 
are truncated at x e . 

As mentioned in fJH resonant interactions are expected 
to dominate the tidal heating. To accurately reproduce these 
resonant interactions, N-b ody simulations need to s atisfy 
several numerical criteria. IWeinberg fc Katzl lj2007aT) pro- 



as that derived for point-mass potentials with vanishingly small 
satellite masses to avoid any further confusion in terminology. For 
extended gravitational potentials, the numerical value will differ 
from 2/3. 



Figure 3. The evolution of fractional mass with time for satellites 
on circular orbits. Curves show the mass fractions remaining for 
satellites truncated at the Roche limit (red), x e (green), and rt 
(blue). All three simulations exhibit continuous mass loss. The 
satellite truncated at the Roche limit loses the most mass and 
the satellite truncated at rt loses the least mass. 



posed explicit requirements for these criteria. First, a suf- 
ficient number of particles are required to cover the phase 
space near resonance (hereafter, the coverage criterion). Sec- 
ond, a sufficient number of particles are required to reduce 
artificial diffusion. Artificial diffusion can come from both 
the gravitational forces of individual particles (hereafter, the 
small-scale noise criterion) and the potential fluctuations 
caused by Poisson noise (hereafter, the large-scale noise cri- 
terion). Besides these particle number criteria, the poten- 
tial solver must also be able to resolve the scale of the res- 
onant potential and the realised phase space distribution 
must cover this region. In our study, we will verify that our 
simulations satisfy all of these criteria. 

Although a simulation can correctly reproduce resonant 
interactions, it is hard to provide a detailed accounting of 
the individual resonances. For intuitive guidance, the reso- 
nant interaction effects can be investigated using perturba- 
tion theor y. We use a numerical pert urbation theory calcula- 
tion as in 



Weinberg fc Katzl (|2007aT ) to investigate resonant 



interaction effects. In this approach, one begins with the nu- 
merical integration of the perturbed orbit-averaged Hamil- 
ton equations in a fixed potential for the entire phase space. 
This step may be followed by an update of the gravitational 
potential. Since this perturbation calculation uses the same 
satellite halo realisation, comparison with the N-body sim- 
ulation result is straightforward. A comparison between the 
results of the N-body simulation and those of the numerical 
perturbation theory calculation provides us with strong ev- 
idence for the existence of resonant dynamics and a definite 
understanding of its role in satellite evolution. 



3 SATELLITE DISRUPTION ON A CIRCULAR 
ORBIT 

To study the resonant torque effect in isolation, we simulate 
a satellite on a circular orbit, which eliminates the gravita- 
tional shock. However, as long as the satellite's rotation does 
not match its orbital rotation, i.e. it is not tidally locked, the 



© 0000 RAS, MNRAS 000, 000-000 



The dynamics of satellite disruption in cold dark matter haloes 5 



satellite will experience a time-dependent azimuthal force 
variation that can result in a net resonant torque. 

3.1 Investigation of the resonance torque effect 

As described in Fig. [3] satellites on a circular orbits contin- 
uously lose mass and have lost a significant fraction of their 
original mass at the end of the simulation (approximately 
the age of the Universe). If a transient adjustment of the 
satellite halo to its host halo potential were the cause of the 
mass loss, the mass loss should lessen and stop in several 
dynamical times as the satellite adjusts to the halo poten- 
tial. Therefore, continuous mass loss from truncated satel- 
lite haloes must result from satellite heating and not from a 
transient initial adjustment. In a circular orbit simulation, 
the magnitude of the external potential does not vary and, 
therefore, there is no gravitational shock. The only mecha- 
nism that can disrupt a satellite under these circumstances 
is a resonant torque. 

The dynamics of the resonant torque are similar to the 
dynamics of the bar-halo interaction. Owing to the coupling 
with the time-dependent perturbation, the orbits whose fre- 
quencies are nearly commensurate with the perturber's fre- 
quency are torqued. The commensurability condition is 



(1) 



where fii and Q2 are the radial and azimuthal frequencies of 
a satellite's orbit, fl per t is the perturber's frequency, and £1, 
I2, and £3 are integers. The orbits that satisfy equation {TJ 
are called resonant orbits and define loci in phase space for 
each triple of integers l\ : I2 : I3. The resonant torque heats 
the satellite, reduces the satellite's binding energy, and en- 
hances its mass loss. Although particular resonant orbits 
receive or lose angular momentum, mass loss is not confined 
to the resonant orbits. Since a satellite is a self-consistent 
gravitational system, the work required to apply the torque 
reduces the entire satellite's binding energy and unbinds the 
most weakly bound halo material. Satellite mass loss in gen- 
eral will be discussed in Sj5] 

Perturbation theory calculations provide helpful physi- 
cal insights into the nature of these resonant dynamics. The 
angular momentum transfer by the resonant interactions 
may be approximated b y a second-order, time-dependent 
pertur bation calculation (|Weinberg||2004l ; IWeinberg fc Kat3 
l2007al lbT). These calculations become prohibitively compli- 
cated for real astronomical systems with multiple time 
scales. First, owing to the finite age of the galaxy and the 
time- dependence of the perturbation, the frequency spec- 
trum becomes broader. Secondly, some resonances need a 
longer time period than the galaxy lifetime to converge into 
the time-asymptotic limit and, in the interim, the first-order 
transient features may strongly affect the response. More 
detailed discussions of perturbatio n theory and resonant 
dynamics are presented elsewhere llTreniaine fc Weinberg! 
1 19841 ; IWeinberg||200i IWeinberg fc Katzl l2007al lbir For our 
calculations here, we represent the host halo potential by 
an orbiting perturber around a stationary satellite with the 
distance of the perturber to the satellite equal to the dis- 
tance between the satellite and the host halo centre in the 
simulation. The mass of the perturber is the enclosed mass 
of the host halo inside the satellite's orbit in the simulation, 
and the perturber's frequency is the same as the satellite's 
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Figure 4. The locations of resonances within the satellite in 
energy-angular momentum space. The low-order resonances with 
— 10 < l\ < 10 and 2 = 1,2, 3,4 are shown (see text for details). 



orbital frequency. We investigate the resonant torque using 
the following three steps: (1) we find the resonances that 
are located within the satellite, (2) we estimate the particle 
number requirements for a given resonance using the crite- 
ria from lWeinberg fc Katzl (|2007al ) ■ and (3) we compute the 
individual resonant effects using a numerical perturbation 
theory calculation and compare them to the simulations. 

Fig. [4] shows the locations of the resonances within the 
satellite phase space. The phase space location is represented 
by energy (E) and normalised angular momentum^ 2 ] (re). The 
range of resonances we examine is — 10 < Zi < 10 and we 
restrict ourselves to / = 1,2,3,4, where I is the spherical 
harmonic of the perturber, in this case the primary galaxy. 
Owing to symmetry, I2 and \h\ should be equal to or less 
than I. Since I2 and I3 are indices of the azimuthal ex- 
pansion in spherical harmonics for the halo and perturber 
respectively, h and £3 shou ld have the same parity as I 
IITremaine fc Weinberg|[l98l ). As one can see in Fig. |4j nine 
resonances are located within the satellite. 

The required numbers of equal mass particles within 
our truncated satellite for the four strongest resonances are 
presented in Fig. [5] There are three particle number criteria 
for each re sonance: coverage, small- scale noise, and large- 
scale noise (|Weinberg fc KatzH2007ah . According to this es- 
timate, we require more than 10 5 satellite particles, which 
translates to more than 3 x 10 5 within the initial virial ra- 
dius, to correctly reproduce the -1:2:2 resonance. This res- 
onance requires the largest number of particles and is also 
the most important resonance as we show below. Hence, 
our 10 6 equal -mass particle halo simulation easily satisfies 
the necessary criteria (see Fig. for all the important res- 
onances. The small-scale noise criterion is not relevant for 
our expansion-code simulations. However, if the N-body sim- 
ulations were to suffer from small-scale noise, such as the 
case in N-body simulations using direct-summation, trees, 
or meshes, an order of magnitude larger particle number 
would be required for these resonances to be modelled cor- 
rectly. 



2 The quantity re = J/ Jmax(E) where J is a particle's angular 
momentum and J m ax(E) is the maximum angular momentum at 
a given energy. 
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Figure 5. Particle number requirements for the -1:2:2, -1:3:3, 
0:1:3, and 1:0:4 resonances. These are the strongest four res- 
onances for a satellite on a circular orbit. Each panel shows 
the thr ee particle number requirements from ["Weinberg Sz Katj 
(2007a) for equal mass particles within the initial truncated satel- 
lite. They are plotted versus r c , the radius of a circular orbit of 
a given energy in units of the host's virial radius. 




Figure 6. The angular momentum transferred from the per- 
turber to the satellite as function of time. The panels show the 
I = 2,3,4 contributions as labelled. The -1:2:2 resonance (top 
panel) dominates the torque. 



Fig. [6] shows the amount of angular momentum de- 
posited by an orbiting perturber calculated using numeri- 
cal perturbation theory for different resonances. The most 
significant angular momentum change is mediated by the 
-1:2:2 resonance. The -1:3:3 and 0:1:3 resonances are the 
strongest among the I3 = 3 resonances, and the 1:0:4 res- 
onance is the strongest resonance among the I3 = 4 res- 
onances. The amount of angular momentum deposited by 
the -1:2:2 resonance is more than an order of magnitude 
larger than the amount of angular momentum deposited by 
any of the Z3 = 3 resonances and three orders of magnitude 
larger than the amount of angular momentum deposited by 
any of the I3 — 4 resonances. Hence, the -1:2:2 resonance 
dominates the resonant satellite torque. 

Fig. [7] shows the location of the angular momentum 
transferred through resonant interactions in phase space by 
plotting the distribution of the change in L z in phase space, 
AL Z , at different times from both the N-body simulation 
and the numerical perturbation theory calculation. We nor- 
malise this change by the total angular momentum at each 



energy E since L(E) at fixed k increases with energy, which 
biases the absolute AL Z distribution toward higher energies. 
By normalising to the total angular momentum at a given 
energy bin, the relative AL Z distribution enhances AL 2 fea- 
tures at low energy. Each panel shows the relative AL Z over 
a fixed time span. Both the phase space location and magni- 
tude of the angular momentum change from the simulation 
agrees well with the numerical perturbation calculation. The 
perturbation theory only includes the resonant torque and 
therefore this agreement is strong evidence that torque by 
resonant interactions is the major mechanism responsible for 
circular orbit satellite evolution. 

Nonetheless, some minor discrepancies remain. First, 
the amplitude of the relative AL Z from the perturbation cal- 
culation is larger than that from the simulation at T =0.0- 
1.0. This results from the abrupt introduction of the pertur- 
bation. In our calculations, we abruptly introduce the exter- 
nal potential (for the simulation) and perturbations by reso- 
nances (for the perturbation calculation) to the initial satel- 
lite. This abrupt introduction may induce a readjustment of 
satellite halo equilibrium and/or cause non-linear features, 
which could cause the discrepancy. Second, the shape of the 
region of angular momentum change is mildly different in 
the simulation compared to the perturbation calculation. In 
particular, the simulation shows a relative AL Z distribution 
in the lower right region of phase space, while the perturba- 
tion calculation does not. The phase space responsible for 
this difference in the AL Z distribution has low binding en- 
ergy and is easily stripped. The simulation results in Fig. [7] 
shows the distribution of the relative AL Z for only the un- 
stripped particles. Unlike in the simulations, the satellite in 
the perturbation calculation does not lose mass. Hence, the 
few minor disagreements between the simulation and the 
perturbation theory are a natural consequence of the ide- 
alisation required to compute the perturbation theory and 
does not invalidate our primary conclusion: resonant torque 
is the major mechanism responsible for satellite disruption 
for a satellite on a circular orbit. 

The absence of a I = 1 resonance contribution is a 
consequence of satellite truncation. Roughly, the tidal ra- 
dius of the satellite coincides with the corotation radius of 
the satellite since Q per t ~ ^2 at r t . The satellite halo is 
isotropic so | °^'' f | < 1 and it follows from equation {1} 
that \h^+l 2 \ < 1 for a resonance within the satellite. Since 
Q1/Q2 is always larger than 1, l^ij^+M has its lowest value 
when Zi = 1 and I2 = —1, making the location of the 1:-1:1 
resonance outside of the satellite. Hence, if the satellite is 
tidally truncated there is no I — 1 resonance contribution to 
satellite disruption. 



3.2 Verifying mass loss by the resonant torque 

The simulation results and the perturbation theory calcu- 
lation agree on the amplitude and phase space location of 
the angular momentum exchange, which confirms the im- 
portance of the resonant torque. However, there is another 
process that may cause mass loss: re-equilibration of the 
satellite after some particles are tidally stripped by the host 
halo potential. This readjustment in the profile forces some 
particles to move beyond the tidal radius and also results in 
continuous mass loss. In our simulation, the satellite initial 
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Figure 7. The distribution of the relative change in L z in phase space for the N-body simulation (left panels) and the numerical 
perturbation theory calculation (right panels) The perturbation theory calculation includes all the resonances in Fig. [4] Except for 
T = 0.0 — 1.0, the results show good general agreement. This agreement provides robust evidence for resonant effects in the circular 
orbit simulation. 
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Figure 8. The mass loss histories from the synthetic force exper- 
iments (see text). 



conditions are designed to reduce the mass lost through this 
process by truncating satellites at the tidal radius but, ow- 
ing to the Eddington inversion, a small amount of mass loss 
is inevitable. To elucidate the affect of resonant torque on 
satellite mass loss, we conduct two types of idealised numer- 
ical experiments. In the first, we add angular momentum to 
the simulation using the torque determined from perturba- 
tion theory, tidal truncation, and the centrifugal and Coriolis 
forces. We call these tests the synthetic force experiments. 
In the second, we simulate satellites corotating with their 
orbital revolution, i.e. as if there were tidally locked. We 
call these the rotating satellite experiments. 

The synthetic force experiments estimate satellite mass 
loss by adding known mass loss mechanisms to a non- 
orbiting satellite simulated in isolation. Satellite particles 
outside the tidal radius do not feel the satellite's gravity, 
allowing particles outside this radius to escape easily. As 
discussed in ij2] we take x e to be the tidal radius, not r t . 
Satellite mass loss can drive structural evolution and change 
$sat(r), hence we update x e during the course of the simu- 
lations. 

We present the results of our synthetic force experi- 
ments in Fig. [5] The curve labelled Tidal (long-dashed) in 
FigHJonly includes tidal truncation and assumes a = in 
the tidal radius calculation (see fj2] and Appendix [S} . We 
include the centrifugal force effect by including it when cal- 
culating x e and we now set a = 0.75 0- We include the 
Coriolis force by adding it to the force used to update each 
particle's velocity at every time step. The mass evolution 
of the resulting simulation is labelled Tidal+ Rotation in 
Fig[8] (short-dashed) . Finally, we include resonant torques 
by adding angular momentum at the rate calculated from 
the numerical perturbation theory calculation plotted in 
FigEl between T = 1.0 and T = 2.0 according to each par- 
ticle's energy and angular momentum. The curve labelled 
Tidal+Rotation+Resonant torque (dotted) includes all three 
mechanisms: tidal truncation, rotation (both the centrifugal 
and Coriolis effects), and the resonant torque. 

We also plot the mass loss from the circular-orbit N- 
body simulation in Fig. [8] (solid line). The idealised simula- 



Note that a is a free parameter and can be between to 1. 



Figure 9. The mass loss history of circular orbit simulations for 
both a non-rotating and a rotating satellite. Since the rotating 
satellite simulation suppresses the resonant torque, its mass loss 
is significantly reduced. 



tion including all the processes matches the actual N-body 
simulation remarkably well. However, the choice of a = 0.75, 
although quite reasonable, was rather arbitrary and we made 
this choice to match the circular orbit simulation. But Fig. 
[8] clearly shows that the resonant torque causes more than 
70% of the satellite mass loss. The remaining 30% of the 
mass loss results from rotational effects. Tidal truncation 
does not cause much satellite mass loss. 

A satellite whose frame corotates with its revolution on 
a circular orbit should not experience a resonant torque. So 
as a further test of the resonant torque mechanism, we com- 
pare the mass loss from a non-rotating satellite simulation 
to a rotating satellite simulation. The non-rotating satellite 
simulation is the same as in H3.ll For the rotating satellite 
simulation, the satellite particles are rotated every time step 
around an axis perpendicular to the orbital plane, which 
passes through the satellite's centre. They are rotated with 
the same frequency as the circular orbit frequency of the 
satellite in the host halo. This fixes the direction of the host 
halo centre in the satellite frame. Owing to this rotation, the 
halo does not feel a time-dependent external potential, and 
hence feels no resonant torque. The two simulations have the 
same initial conditions. However, owing to the artificial rota- 
tion, the rotating satellite feels an extra Coriolis force. Dur- 
ing the course of rotating satellite simulation, we subtract 
this extra force to make the comparison more meaningful. 

Fig-EUshows that the rotating satellite loses only » 40% 
of the mass lost by the non-rotating satellite by T = 5.0. 
Moreover, the mass loss histories of both satellites are almost 
identical at early times (T < 0.7). This mass loss results 
from the adjustment of the satellite's initial conditions to the 
host halo potential. If we ignore satellite mass loss during 
this phase, the non-rotating satellite loses about five times 
more mass than the rotating satellite, again confirming that 
resonant torques are the dominant mechanism responsible 
for circular orbit satellite mass loss. 

The mass loss rates in the Tidal+ Rotation simulation 
in Fig.|S]and the Rotating satellite in Fig.|5]are almost iden- 
tical if we ignore the early relaxation phase in the rotating 
satellite simulation. Both simulations attempt to reproduce 
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Figure 10. The circular, eccentric, and inner orbits. Top panel: 
the three trajectories in the orbital plane. Bottom panel: the time 
evolution of the distance between the satellite and the host halo 
centre. 

circular orbit satellite mass loss without resonant torque ef- 
fects, using totally different approaches. This similarity con- 
firms that our treatment for non-resonant mass loss captures 
the essential dynamics and implies that that our estimates 
of satellite mass loss in these experiments is reasonable. In 
conclusion, Figs.|S]and|5]confirm that resonant torque effects 
dominate circular orbit satellite mass loss. 



4 SATELLITE DISRUPTION ON AN 
ECCENTRIC ORBIT 

Several physical processes, including resonant torques, grav- 
itational shocks with resonant shock effects, and continuous 
tidal truncation, simultaneously affect a satellite on an ec- 
centric orbit. To study this complicated interplay and to un- 
derstand the characteristics of the individual physical pro- 
cesses, we compare the evolution of three different simula- 
tions: a circular orbit simulation, an eccentric orbit simula- 
tion, and an inner orbit simulation (see Fig. [TDJ. The cir- 
cular orbit simulation is that described in ij3] The eccentric 
orbit simulation is an e = 0.5 orbit with an apocentre of 
r ap o = 0.6 and a pericentre of r per i = 0.2 Q. The circular 
orbit and the eccentric orbit have the same orbital energy. 
The inner orbit simulation is an e = 0.72 orbit with an 
apocentre of r apo = 0.4 and a pericentre of r per i = 0.064. 
The normalised angular momentum (k) of the inner orbit is 
0.55, which is near the median k of subhaloe s in a sample 
taken from recent cosmo logical simulations (|Ghigna et al.l 
ll99Sl : IZentner et al.lliooBh . 

Fig. [11] shows the evolution of the enclosed mass within 
90%, 22%, 2.2%, and 0.22% of the original satellite radius. 
The two eccentric orbit satellites lose their mass at peri- 
centre episodically, while that for the circular orbit appears 
continuous. The eccentric-orbit satellites lose significantly 
more mass than the circular orbit satellite. 

The computation of a satellite's gravitationally bound 
mass is difficult in practise and we use the mass enclosed 

4 We define eccentricity as e = (r apo — r peri ) / (r apo + r peri ) 
© 0000 RAS, MNRAS 000, 000-000 



Figure 11. The mass evolution for the circular, eccentric, and 
inner orbit simulations. We plot the mass within radii of 0.045, 
0.04, 0.01, and 0.001 versus time in the top left, top right, bottom 
left, and bottom right panels, respectively. The initial satellite 
radius is 0.045. 
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Figure 12. Evolution for the circular, eccentric and inner or- 
bit simulations. Top panel: the mass enclosed within the initial 
satellite radius. Middle panel: the mass enclosed within the in- 
stantaneous tidal radius. Bottom panel: the evolution of the tidal 
radius. 



within the truncation radius as a proxy. Fig. [12] shows the 
mass loss history for the three simulations using two different 
definitions of the bound satellite mass. The first definition is 
the total mass inside the initial satellite radius. The second 
definition is total mass inside the instantaneous tidal radius. 
We find that a satellite's mass rapidly decreases and then 
increases again around the time of pericentre passage, if we 
use the tidal radius definition for the bound mass. This mass 
evolution results from the rapid variation in the instanta- 
neous tidal radius (Fig. 1121 bottom) . However, if this radius 
varies faster than the orbital time of an escaping particle, 
the tidal radius will not represent the dynamics of escape. 
Therefore, we generally use the initial satellite radius defini- 
tion for our working definition of the bound satellite mass. 
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Figure 14. The figures show the enclosed mass evolution of two 
simulations: the fiducial simulation, which is the same as the inner 
orbit simulation, and the 'turn-off' host simulation, where the 
host halo potential is turned off at pericentre. 



4.1 Heating by the gravitational shock at 
pericentre 

According to Fig. 1111 satellites on eccentric orbits lose mass 
mostly at pericentre. This suggests that satellite heating by 
the gravitational shock at pericentre is the dominant mecha- 
nism for satellite mass loss and disruption. We can discrimi- 
nate two distinct effects of this gravitational shock: evolution 
of the satellite's structure and the persisting tidal trunca- 
tion. 

The top panel in Fig. H3l shows overlapping snapshots of 
the inner orbit simulation over one orbital period beginning 
from apocentre. The colour scale is logarithmic in satellite 
density. The lower-right panel shows the satellite evolution 
during pericentre passage. The satellite is compressed per- 
pendicular to its trajectory and stretched out along its tra- 
jectory by the external potential. This deformation causes 
the density to increase, as observed in Fig. [11] where the 
enclosed mass within 90% (the top right panel), 22% (the 
bottom left panel), and 2.2% (the bottom right panel) of 
the original satellite radius increases at pericentre. The two 
bottom panels in Fig. [13] show blown-ups of the two boxed 
areas in the top panel. To highlight density differences, the 
colour now scales linearly with satellite density from p = 1 
to p = 100, and white represents p > 100. The area of the 
white region in the right panel is larger than that of the left 
panel because of the density enhancement at pericentre. 

The energy input from the gravitational shock reduces 
the satellite's binding energy, causing it to expand, and leads 
to tidal stripping. Work is done on the satellite potential 
to strip particles, this extra energy drives expansion and 
enhances the mass loss. In addition to enhancing the mass 
loss, the tidal truncation causes the satellite haloes to lose 
equilibrium. Fig. [12] shows that when the satellite moves 
deep inside the host halo, the tidal radius decreases. A large 
fraction of the satellite material outside of the tidal radius is 
not stripped and reenters the satellite when it again moves 
outward in the host halo and the tidal radius again increases. 
This rapidly changing tidal radius keeps the satellite haloes 
out of equilibrium. 

To help understand the effect of the host halo potential, 



we compare two inner orbit simulations. The first simulation 
(fiducial) is the same as the original inner orbit simulation. 
The second simulation (turn-off host) is the same inner or- 
bit simulation but we artificially "turn off" the host halo 
potential at pericentre (T w 0.45). We present the mass 
loss histories of these two simulations in Fig. 1141 After turn- 
ing off the host halo potential, the mass loss history of the 
two satellites diverges. The turn-off host simulation shows 
~ 15% less mass loss than fiducial simulation. In addition, 
the satellite mass loss history of the turn-off host simula- 
tion shows an increase in mass after the significant mass 
loss caused by the gravitational shock, while the mass loss 
history of the fiducial simulation shows a continuous mass 
loss after the shock. The increase in mass results from lost 
particles returning to the satellite. Conversely, the persist- 
ing tidal truncation by the host halo potential in the fiducial 
simulation enhances satellite mass loss, causing the satellite 
halo to lose equilibrium and makes the evolution nonlinear. 

In summary, both the internal structural evolution and 
the subsequent tidal truncation driven by time-dependent 
heating leads to significant evolution. These nonlinear phys- 
ical processes are essential ingredients for an accurate pre- 
diction of satellite disruption and, therefore, simple approx- 
imati ons such as impulsive heating are insufficient |Spitzerl 
Il987t) . 



4.2 Resonant heating on an eccentric orbit 
satellite 

We have divided the resonant dynamics of satellite evolution 
into a shock and a torque, but this distinction is arbitrary in 
principle. For example, both resonant effects simultaneously 
drive the evolution of a satellite on an eccentric orbit. A 
consistent perturbation theory calculation for an eccentric 
orbit perturber incorporates the full time dependence, cou- 
pling both the radial and the azimuthal orbital frequencies 
to a satellite's phase space. In addition, orbital decay adds 
the complication of a continuous spectrum of perturbation 
frequencies. Hence, a perturbation theory calculation for an 
eccentric orbit is a difficult and expensive. To get around 
this difficulty, we estimate the response by considering the 
most important frequencies one at a time. The strongest 
resonant effects are associated with the pericentric passage. 
We can, therefore, approximate the coupling with that of a 
circular satellite orbit at the pericentre radius and proceed 
as in i|3.1l 

Fig. ll5l shows the locations of the resonances in satellite 
phase space for this equivalent circular orbit problem. We 
examined resonances with h £ [—10,10] and I £ [1,4]. In 
addition to the -1:2:2 resonance, the 0:1:1, 1:-1:1, and 1:0:2 
resonances occur within the satellite's phase space. Com- 
paring this to the results from the larger radius circular 
orbit calculation f H3.ip . the coupling of the larger instan- 
taneous angular frequency at pericentre to higher binding 
energy orbits that have larger orbital frequencies, moves the 
resonances inward in radius. Hence, some resonances that 
were previously outside the satellite in i|3.1l are now within 
the satellite. In addition, the stronger tidal field at smaller 
orbital radii in the host increases the strength of the cou- 
pling. Using numerical perturbation theory, the strength of 
the 1:-1:1 and the 1:0:2 resonances are now comparable to 
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Figure 13. Snapshots of the inner orbit simulation for one radial period beginning at apocentre. Top panel: Overlapping snapshots of 
the satellite, shown at T =0.0, 0.16, 0.31, 0.47, 0.63, and 0.78. The origin of the inner axes (green) is the centre of the host halo. The 
dark matter density scales logarithmically with colour. Bottom panels: The satellites at T = 0.0 and T = 0.31 zoomed in to the size of 
the small boxes (green) in the top panel. To emphasise the dark matter density differences, the colour scale is now linear in density. The 
bottom snapshots show an increase in the size of the high density regions at T = 0.31 at pericentre compared to the satellite at T = 0.0 
at apocentre. 
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Figure 15. As in Fig. [4] but assuming a satellite frequency of 
£lsat = 6.56, which corresponds to the frequency of a circular 
orbit at r per { for the eccentric orbit simulation. Only the 1 = 1 
and 1 = 2 resonances are shown. Unlike in Fig. [4] the 1:-1:1, 0:1:1, 
and 1:0:2 resonances are located inside the satellite. 



the strength of the -1:2:2 resonance, which is much stronger 
than the -1:2:2 resonance in the larger radius circular orbit. 

In summary, resonant heating is enhanced for an eccen- 
tric orbit both because of the larger number of resonances 
within the satellite and the stronger tidal force felt by the 
satellite at pericentre. However, our approximate perturba- 
tion theory calculation does not include the full time depen- 
dence of the eccentric orbit, a calculation that is extremely 
difficult and is beyond the scope of this paper. 



5 SATELLITE STRIPPING AND ITS EFFECTS 
5.1 The stripping process 

Owing to satellite heating, particles in a satellite gain en- 
ergy and angular momentum. This reduces the satellite's 
binding energy and enhances the satellite's mass loss. There 
are several mechanisms that heat a satellite and each one 
is effective over different ranges of binding energy. However, 
the re-equilibration of the satellite tends to globally redis- 
tribute the work throughout the satellite profile, washing 
out the nature of its origin. The observational signatures of 
the subsequent mass loss, then, tend to be universal. 

Fig. [16] plots the fraction of bound particles remain- 
ing in different parts of phase space at different times for 
a satellite on a circular orbit. Low binding energy particles 
are stripped at earlier times and high binding energy parti- 
cles are stripped later. In other words, the satellite stripping 
process is an outside-in process in energy space. We can un- 
derstand this behaviour as follows. Weakly bound satellite 
orbits are affected and stripped by the tidal force beyond a 
characteristic radius. The apocentres of orbits with a given 
energy are within a factor of two of the radius of a circular 
orbit with the same energy, even for a zero-angular momen- 
tum orbit. For example, Fig, ll7l marks the radius of a circular 
orbit (ft = 1) and the apocentre of a radial orbit (ft = 0) 
for an energy of E — —0.044; the apocentre of the radial or- 
bit is only 50% larger than the radius of the circular orbit. 
Therefore, energy and not the relative angular momentum ft 



Figure 16. The fraction of particles remaining in phase space 
at different times. The colour code is linear from all particles 
stripped (0) to no particles stripped (1). Each panel shows the 
fraction of particles remaining at the labelled time compared to 
the initial satellite. Satellite stripping is an outside-in process in 
energy space. 
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Figure 17. The apocentre and pericentre of orbits in the initial 
satellite model. Dashed lines represent the apocentre and solid 
lines the pericentre for orbits with energies of E =-0.076, -0.044, 
and -0.027. The radii for circular orbits with these energies are 
0.01, 0.02, and 0.03, respectively. The x-axis is the distance from 
the halo centre and the y-axis is the normalised angular momen- 
tum, k. The area between the pericentre and apocentre lines is 
the range in radii for orbits with each E and k. The two arrows 
show the circular orbit radius and the apocentre radius for a ra- 
dial orbit (k = 0) with E =-0.044. 



determines the stripping boundary. Although very low angu- 
lar momentum orbits can be stripped at lower energies, Fig. 
[TB] shows that any trend towards a larger escape fraction for 
smaller k at fixed energy is very weak. 

We plot the fraction of particles remaining in different 
parts of phase space for the eccentric orbit simulation in 
Fig. 1181 The stripping process remains an outside-in process 
in energy space but now the trend toward a larger escape 
fraction at smaller relative angular momentum is stronger. 
This angular momentum dependence results from the time 
dependence of the tidal radius as shown in Fig. 1121 Near 
a satellite's pericentre, satellite particles with low angular 
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Figure 18. The same as Fig. [16] but for the eccentric orbit sim- Figure 19. The circular velocity profile of the LMC dark matter 

ulation. halo at the initial and present time. 



momentum are temporarily outside the instantaneous tidal 
radius and are prone to escape. 

The outside-in stripping in energy space has interesting 
consequences for satellite evolution. Since the stars and cold 
gas in a satellite galaxy are centrally concentrated in their 
own dark matter halo, it protects them from stripping. The 
degree of protection for any component depends on its bind- 
ing energy. For example, consider a disk and elliptical galaxy 
with the same linear exten10 and the same orbit in the host 
halo. The disk galaxy contains material mostly on circular 
orbits, which have the largest binding energy at fixed radius, 
while the elliptical galaxy is made from a wide range of stel- 
lar orbits. We therefore predict that stars in the elliptical 
galaxy will be stripped sooner than stars in the disk galaxy. 
This difference should be incorporated into models of satel- 
lite galaxy stripping and when one estimates the subsequent 
stellar halo distribution. 



5.2 The LMC stellar tail 

The LMC is the largest Milky Way satellite. It has been the 
subject of numerous observational studies and characteris- 
tics of its structure and kinematics have been well estab- 
lished. However, the origin of the Magellanic Stream, a thin 
neutral hydrogen tail stretching over 100° along a Galactic 
great circle, remains unexplained. Observations suggest that 
the Magellanic Strea m may be a relic of a past interaction 
with the Milky Way (jMathewson et alJlraTu ; IPutman et all 
1998). If this is the case, tidally stripped stellar ejecta should 
be observed as predicted by tidal interaction theory. How- 
ever, none of the many searches for stellar ejecta has suc- 
cessfully detected a population of stars connected with the 
Magellanic Stream. Several scenarios have been suggested to 
resolve this conflict. First, the LMC gas may be removed by 
ram- pressure forces caused by an interaction with the outer 
disk llMoore fc Davislll994ri or with hot gas in the Galactic 
halo ( Mastropietro et all 20051 ) instead of a tidal interaction. 
Second, the Stream may be a remnant from the tidal dis- 
ruption of the Small Magellanic Cloud (SMC) owing to its 



5 In fact elliptical galaxies are generally more concentrated than 
disk galaxies. 



gravitational interaction with the LMC and t he Milky Way 
|Murai fc Fuiimotolll98fj| ; IConnors et aLlbOOel ). Based on the 
results from the previous sections, we suggest that the LMC 
stars are protected by its dark matter halo while the Stream 
is the remnant of a low binding energy, extended neutral gas 
disk, analogous to the Milky Way's extended HI disk. Al- 
though smaller satellites and globular clusters have symmet- 
ric tidal tails, which has been used as an a rgument against 
a tid al origin for the Magellanic Stream (|Moore fc Davisl 
1994), larger satellites have a asym metric tidal tails owing 
to the asymmetry in the tidal force IjChoi et al.1120071 ) . After 
the gas is stripped it would then hydrodynamically interact 
with th e hot halo gas, changing the orbit of the gas in the 
stream l|Mastropietro et al.ll2005T ). 

We implement an idealised simulation of LMC evolu- 
tion in the Milky Way dark matter halo to test this hy- 
pothesis. The static Milky Way dark matter hal o potential 
is based on the Al model of iKlypin et al. (2002); a c = 12 
NFW halo with R vir = 258 kpc, M vir = 1.0 x 1O 12 M , 
and V m ax = 163 km/s. The live LMC halo model, which is 
also based on a c = 12 NFW halo model, is truncated at 
its x e radius assuming that the LMC is located at 0.6R V i r . 
To make a stable halo model we performed an Eddington 
inversion to the truncated LMC halo model, as presented 
in Sj2] and we realise a phase space of 10 6 equal-mass par- 
ticles. The initial LMC halo circular velocity, Vmax,LMC = 
QA2V max ,host, is chosen to match the obse rved LMC circu- 
lar velocity of V ma x lmc = 0.3V max ,host l|Kim et al.l Il998l : 
Ivan der Marel et alj|2002h at the present day (see Fig. I19[) . 

Since the LMC's evolutionary history is not well- 
defined, we assume a simple LMC orbit: an e = 0.5 
orbit with a pericentre at 0.2i?„; r and an apocentre at 
0.6-R^i,-. The curr ent LMC is assumed to be located just af- 
ter it s pericentre l|Gardiner et al.| [T994; Gar diner fc Noguchil 
1996). The simulation starts at 0.6R V i r and runs a little 
more than two radial periods to represent the present-day 
LMC. It is worth noting that a consensus on the LMC or- 
bital parameters has not yet been reached. Recent proper 
motion measurements u sing the Hubble Space Telescope 
l|Kallivavalil et all l200fj ) imply a very large galactocentric 
velocity. This large LMC velocity combined with an im- 
proved Milky Way mass model suggests either that the LMC 
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Figure 20. The distribution of particles remaining in the LMC 
halo at different times. Left panels: The fraction of particles re- 
maining in the LMC halo based on their initial location in phase 
space (see Fig, ll6l l, Bight panels: The mean radius of the particles 
remaining as a function of phase space coordinates. 



5.3 Density profiles and rotation curves 

The left panels in Fig. [21] show the evolution of the den- 
sity profiles for the three different orbit simulations. As ex- 
pected given the mass loss rates shown in Fig. 1111 the inner 
orbit satellite's density profile evolves the most and the cir- 
cular orbit's density profile evolves the least. However, even 
though the mass loss rates are different, the overall trends in 
the density profile evolution for the three satellites are sim- 
ilar. A decrease in density occurs at all radii, which is par- 
ticularly evident in the inner orbit simulation, where even 
at r < 0.001 the density decreases uniformly with time. Al- 
though the density decreases, the steepness of the density 
profile does not decrease. Hence, the inner cusp is preserved 
during tidal mass loss, i.e. the evolution of satellites owing to 
interactions with their host halo will not remove NFW cen- 
tral cus ps to make them agree more with observe d density 
profiles (jDiemand et al ] |2007l : ISpringel et afl|2008h . 

One can see the effect on the circular velocity profiles for 
the satellites on the three different orbits in the right panels 
of Fig. 1211 Once again, even though the peak circular veloc- 
ity drops from, the shape of the circular velocity profile does 
not change dramatically. In particular , the central c ircular 
velocity profile is still steeply rising. IStoehrl (|2006h claim 
that the circular velocity profiles of subhaloes are best fit by 
a parabolic function, which is shallower than NFW circular 
velocity profiles, and that this confirms their earlier result 
that the internal structure and kinematics of the Milky Way 
satellites are in goo d agreement with t he subhaloes found in 
CDM simulations IjStoehr et al.l 120021 ). The persistence of 
the steep inner density structure in our satellite simulations 
conflicts with their assumed parabolic velocity profile. Re- 
cent very high resolution simulations following the formation 
of CDM haloes in a cosmological context also have steep cen- 
tral de nsity profiles in their subhaloes, which co nfirms our 
result (|Diemand et alj|2008l ; [Springel et alj|2008h . 



is on its first passage about the Milky Way or that its orbital 
period and apocen tre must be a fac tor of 2 larger than pre- 
viously estimated (jBesla et al.ll2007l ). In this paper, we take 
the classical LMC orbital p arameters, which are bas ed on 
the prescription induced by iMural fc Fuiimotol (|l980h . If it 
turns out that the classical point of view is false, our quanti- 
tative estimation should be revisited. However, our qualita- 
tive interpretation about the pattern of the LMC disruption 
is not sensitive to our particular choice of LMC orbital pa- 
rameters. 

Fig. 1201 shows the fraction of particles remaining in the 
LMC in different parts of phase space (left panels) and their 
average distance from the LMC centre (right panels) at dif- 
ferent times. Once again, the stripping corresponds to an 
outside-in process in energy space. Fig. [20] shows that parti- 
cles with E < —0.5 have not yet been stripped. We assume 
that the LMC stars are mostly on circular orbits (k > 0.9), 
because the visible LMC is a rotationally supported sys- 
tem. The current mean radius of the particles for E = —0.5 
and k = 0.9 is about 15 kpc, and hence LMC stars inside 
15 kpc are expected to be protected by the dark matter 
halo, which is larger th an the observed disk radius of 10 
kpc (jKunkel et al.l 119971 ). This suggests that the detectable 
signature of the LMC stellar ejecta could be very small. 



6 ESTIMATING SATELLITE MASS LOSS 

To estimate satellite mass loss, researchers often use the im- 
pulse approximation. The impulse approximation assumes 
that the perturbation time sc ale is much s horter than the 
internal dynamical time scale (Spitzcr 1987). However, in fJ2] 
and 33] we demonstrated that satellite halo evolution results 
from resonant effects, internal structure evolution, and tidal 
truncation in addition to impulses. In this section we will 
include the two resonant effects described in ij3]to produce 
an improved analytic estimate. We implement these approx- 
imations in a simple model that computes the mass loss by 
tidal truncation, gravitational shocks, and resonant torques 
in spherical shells. For the tests presented here, we use the 
same satellite models and orbits presented in §32HH At each 
time step, we proceed as follows: 

(i) The tidal radius x e is computed for the current posi- 
tion in the satellite's orbit. The instantaneous satellite mass 
is the mass of all particles with r < x e and the remainder is 
considered to be stripped. 

(ii) The satellite's mass outside of the tidal radius is grad- 
ually stripped on a crossing time scale. In particular, for each 
time step we remove a mass fraction of At/torb, where t or b 
is the orbital period at the half mass radius of the satellite. 
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(c) The inner orbit 



Figure 21. The density profiles (left panels) and rotation curves (right panels) of the satellite dark matter haloes for the three satellite 
orbit simulations at various times. 



This procedure prevents unrealistic immediate and perma- 
nent satellite mass stripping outside the tidal radius (see 

S3)- 

(iii) The changes in energy and angular momentum for a 
given mass shell caused by the time-dependent external po- 
tential are computed in two parts: 1) the change in energy 
from gravitational shocks is only computed at every peri- 
centre; and 2) the change in angular momentum owing to 
the resonant torque is computed at every step. 

(iv) An increase in energy and angular momentum gener- 



ally results in reducing the density of a mass shell. Overall, 
the reduction in density results in mass loss, decreases x e , 
and hence drives further mass loss. 

This mass loss algorithm is simila r to the scheme used in 
iTavlor fc Babull l|200ll . 12004 l2005al lbh. It provides the satel- 
lite mass loss history, which can be used to compare with 
the simulation results. 

We determine the tidal radius as the satellite radius 
where d$ e // = (see ^2]and Appendix |A")) . Since this model 
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does not include any velocity information about individual 
satellite particles, our estimate cannot include the Coriolis 
force. Instead we choose a — 1.0 when evaluating $ e // to 
compensate for ignoring the Coriolis force. 

For the gravitational shock heating, we compute the 
change in energy using the impulse approximation with an 
adiab atic correction. Ass uming an spherical isothermal host 
halo. lGnedin et all (|l999h derive the following expression for 
the energy change: 



AE 



R%V V 



M halo 
R v 



A(x) 



(2) 



where r is the satellite's radius, R p is the pericentre distance, 
and V p is the satellite's velocity at pericentre. Mhalo is the 
mass of the perturber, which for us is the enclosed halo mass 
at pericentre. The quantity A(x), the adiabatic correction, 
reduces the heating if the internal orbital time is significantly 
shorter than the impulsive time scale. We explor e both the 
Spitze r correction and the Weinb erg correction (|Weinberd 
Il994al lbl; IGnedin fc Ostrikerlll999T l : 

2\ 



A(x) 
A(x) 



exp(—2x 
(1 + x)- 1 - 5 



(Spitzer) (3) 
(Weinberg) (4) 



where x — ujtvR p /V p and ui is the azimuthal frequency 
of an orbit at r. Equations ©-(SI) determine the energy 
change for each mass shell and for each shock event. The 
added energy expands the satellite. Using the virial theo- 
rem, the expansion ra te can be estimated as Ar = AEr 2 
l|Tavlor fc Babulll200ll ). This satellite expansion reduces the 
density in each shell: Ap = -(9/2n)(AE/r 2 ). We use this 
relationship to compute the mass loss in Step (iv) of the 
algorithm above. 

In addition to heating by gravitational shocks, the re- 
sults from ^confirm that heating by resonant torques plays 
an important role in satellite mass loss. This effect has not 
been included in most studies of satellite evolution. To accu- 
rately estimate the reso nant torque, one needs h igh quality 
N-body simulations (see Weinberg & Katz 2007a) or pertur- 
bation theory calculations IjWeinbe rg 1986, 1989). However, 
we can estimate the torque using a simple model as follows. 
First, we compute the torque on a mass shell by assigning 
the mass in the shell to an equivalent ring. The torque on 
the ring is computed using a simple spin-orbit coupling cal- 
culation as in a planet-satellite interaction. In this approx- 
imation, we assume that the host halo potential is a point 
mass at the centre of the host halo. We model the quadru- 
ple moment of a satellite mass shell as two equal, diamet- 
rically opposed point masses in the satellite's orbital plane, 
separated by the diameter of the mass shell. Owing to the 
separation, the gravitational forces on the two masses differ. 
The differential gravitational force torques the two masses 
and we evaluate the total torque on the satellite's mass shell 
by integrating this torque through the ring. 

The mathematical form of the torque on the mass shell 
is then 

r rlng {R) = 6GM-5- (5) 

Rsat 

where R sa t is the distance from the host halo centre to the 
satellite, R is the distance from the satellite centre to a 
given mass shell, and M is the enclosed host halo mass 
(Murray fc DermottJ Il999l . Chapter 5.3 see their Fig. 5.7). 



No Shock 
Torque 
Simulation 



Figure 22. The mass loss history of the circular orbit simulation. 
Analytic estimates of the mass loss history are compared with the 
simulation. No Shock is an analytic estimate only including tidal 
truncation and Torque is an analytic estimate including the effects 
of tidal truncation and resonant torques. 



Equation © provides the specific torque on the satellite 
mass shell. 

Second, we compute the fraction of resonant orbits in 
each mass shell. Equation ([5]) computes the torque on a rigid 
mass shell, not the torque on the mass shell owing to the res- 
onant interaction. To estimate the resonant torque, we esti- 
mate the mass fraction of resonantly coupled orbits in the 
mass shell and multiply the torque on the rigid mass shell 
by this fraction. As shown in equation ^Q, determining res- 
onant orbits requ i res det ai led calculations. Us i ng eq uation 
(24) of I Weinberg! (|l994bh . IGnedin fc Ostrikerl (l999) show 
that the number of stars at the peak amplitude scales as 
1/t, where r is the characteristic duration of the shock. The 
resonant shock and the resonant torque are based on the 
same physics with a different coupling and frequency. By 
dimensional analysis, 1/t ~ uj, where uj is the angular fre- 
quency of satellite mass shell Q. Therefore, we estimate this 
fraction as the ratio of the azimuthal angular frequency of 
the orbiting satellite to the angular frequency of the satel- 
lite's mass shell. Consequently, the torque on the satellite 
mass shell becomes: 



T S at(R) = Tring{R) X < 



u>(R) 



SI, 



if Q s < oj{R) 



if Q s > lj(R) 



This torque leads to expansion and the expansion rate 
be estimated as 



Ar = 



4r 



M sat (r 



■TsatAt. 



(6) 



can 



(7) 



This expansion reduces the satellite's density as Ap oc — Af- 
We also use this relationship to compute the mass loss in 
Step (iv) of the algorithm above. 



6 The angular frequency u> = -S where V c is the circular velocity 
and R is the radius of a given mass shell. 
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No Shock 
Shock(S) 
Shock(W) - 
Shock(W)+Torque 
Simulation 



Figure 23. The mass loss history for the eccentric orbit sim- 
ulation. Analytic estimates of the mass loss history are com- 
pared with the simulation results. No Shock is an analytic es- 
timate only including tidal truncation. Shock(S) is an analytic 
estimate including tidal truncation and impulse shock heating 
with the Spitzer correction. Shock(W) is an analytic estimation 
including the effects of tidal truncation and impulse shock heating 
with the Weinberg correction. Shock(W) + Torque is the same as 
Shock(W) but with the addition of the resonant torque. 



No Shock 
Shock(S) 
Shock(W) - 
Shock(W)+Torque 
Simulation 



Figure 24. The same as Fig. I23l but for the inner orbit simula- 
tion. 



rithm to the simulations. Here, we include the No Shock 
model, the impulse approximation with the Spitzer correc- 
tion [Shock(S)], the impulse approximation with the Wein- 
berg correction [Shock(W)], and the impulse approximation 
with Weinberg correction together with our resonant torque 
approximation [Shock(W)+ Torque]. All the estimates in- 
clude tidal truncation. The Shock(W) + Torque model best 
represents the mass loss seen in the simulation while the 
Shock(W) with no resonant torque predicts less mass loss 
than that seen. The importance of the resonant torque is 
most obvious for the inner orbit simulation. Compared to 
the simulation, the estimate including the shock and the res- 
onant torque [Shock(W) + Torque] does significantly better 
estimating the mass loss in the simulation than that of the 
shock alone, which significantly underpredicts the mass loss. 
Tidal truncation alone [No Shock] shows the worst agree- 
ment. 

Although our new algorithm provides a dramatically 
improved satellite mass loss history, discrepancies remain 
owing to the many complicated physical processes involved 
in satellite disruption. First, the compressive gravitational 
shock is a strongly non-linear, time-dependent perturbation. 
Although the adiabatic correction includes some of this time 
dependence by including the work done at the frequency 
peak, many channels of possible coupling are ignored. Sec- 
ond, satellite heating by resonant interactions is too compli- 
cated to be accurately represented by our simple parametri- 
sation. For example, the adiabatic correction extends the 
coupling to lower energies by noting that the resonant cou- 
pling will have a power law rather than an exponential scal- 
ing with frequency. However, for a particular interaction, 
individual resonances may dominate the response; this is 
analogous to the difference between a line and continuous 
spectrum. Third, a real satellite halo has a wide range of 
orbits, from circular to radial. The mass shell scheme can- 
not accurately capture the dynamics of these different or- 
bits. Moreover, the entire satellite structure readjusts its 
evolution in the course of satellite disruption. This readjust- 
ment process is also not included. These arguments suggest 
that although our new algorithm provides improved satellite 
mass loss histories, high quality simulations are still neces- 
sary for an accurate prediction of mass loss. 



The mass loss histories for the circular orbit simulation, 
the eccentric orbit simulation, and the inner orbit simulation 
are compared with the results of our mass-loss algorithm 
in Figs. [22] - 1241 respectively. For the circular orbit case, 
there are no gravitational shocks. Fig. !22l shows two analytic 
estimates: No Shock and Torque. The No Shock estimate 
only includes tidal truncation and predicts negligible mass 
loss. The Torque estimate includes both the resonant torque 
approximation and tidal truncation. Its mass loss history 
is much more similar to the simulation: more than 50% of 
the original mass is lost. This suggests that mass loss for 
a satellite on a circular orbit mainly results from resonant 
torques and our algorithm for resonant torque provides a 
dramatically improved description. 

For the eccentric orbit and inner orbit simulation cases, 
the gravitational shock plays an important role in driving 
mass loss. Figs. [23] and [24] compares our mass loss algo- 



7 SUMMARY AND CONCLUSIONS 

Using high resolution simulations with cosmologically mo- 
tivated initial conditions, we investigate the physical pro- 
cesses responsible for the evolution of satellite galaxies in 
their host halo. We identified and explored the following im- 
portant physical mechanisms that result in satellite galaxy 
disruption. Our main results are as follows: 

(i) Resonant mechanisms of two types play a key role 
in satellite disruption: the resonant shock and the resonant 
torque. 

(ii) We studied satellites on circular orbits to isolate the 
effects of resonant torques. The ILR-like resonance l\ : h : 
^3 = — I : 2 : 2 (see equation |T}) dominates the torque 
resonance because all the I = 1 resonances and the 0:2:2 
(corotation-like) and the 1:2:2 (OLR-like) resonances are lo- 
cated outside the satellite's tidal radius. 
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(iii) Some important resonances require more particles 
than are typically used in cosmological simulations. For ex- 
ample, to accurately reproduce the heating from the -1:2:2 
resonance, a satellite simulation needs more than 10 par- 
ticles within the virial radius for our expansion code and 
possibly more for codes affected by small-scale noise, e.g. a 
tree code. Too few particles results in less mass loss and it 
suggests that the lifetime of dark matter subhaloes in cur- 
rent cosmological simulations could be overestimated. 

(iv) For satellites on eccentric orbits, gravitational shocks 
dominate the heating but the heating rate is significantly 
underestimated by the impulse approximation for several 
reasons. First, strong low-order resonances couple to the 
phase space in the 'adiabatic' regime and the work done 
on the satellite drives structural evolution. Second, the on- 
going tidal truncation affects the dynamics of the escaping 
material and the satellite's size. The interplay of these mech- 
anisms leads to continuous mass loss and evolution. 

(v) The number of resonances in a satellite increases as 
the satellite's orbital frequency increases; the strength of 
each resonant interaction also increases. This affects satel- 
lites on low-energy and low-angular momentum orbits which 
have small pericentres and high satellite azimuthal frequen- 
cies. 

(vi) Satellite mass stripping is an outside-in process in 
energy space. We discuss the morphological consequences 
for the stripping signatures of stellar and gaseous streams. 

(vii) We present an improved algorithm to estimate satel- 
lite mass loss including both resonant shocks and torques. 
We show that the mass-loss history computed using this 
algorithm reproduces the general features seen in the simu- 
lations. 



The satellite 
has b een studied 
fe.g. | Havashi et al 



evolution using N-body simulation 
by a n u mber of authors re cently 
>003) ; iKazantzidis et al.l 12004 ; 
iBovlan-Kolchin fc Mai 120071 ). Their approach differs 
ours in several ways. We use a highly idealised simulation 
configuration to better discern the details of the dynamical 
mechanisms, focussing on the detailed physical processes 
affecting satellite evolution such as resonant dynamics that 
have not been rigorously addressed elsewhere. For example, 
we equilibrated the initially tidally truncated satellite to 
reduce the artificial perturbation from sudden introduction 
of the host halo potential. We also explored a circular orbit 
simulation to clearly demonstrate the effects of resonant 
torque. We performed a linear perturbation calculation to 
estimate the required resolution for resonant effect, and 
made the simulations satisfy this requirement. In addition, 
we use an expansion code to reduce the force fluctuations 
on small scales. Owing to these efforts, we were able to 
demonstrate the importance of the resonant torque in 

sate llite disruption . 

lHavashi et all (|2003T ) used an NFW halo as a satellite 
initial conditions and a Tree code for a potential solver. They 
found that the simple tidal-limit approximation underesti- 
mates the mass loss, as we do, an d found structural evolu - 
tion , also similar to our findings . IKazantzidis et aT] (|2004T ) 
and lBovlan-Kolchin fc Mai l|2007h also found that similar in- 
ternal structure evolution using tree-code simulations. This 
consensus suggests that the inner cusp of satellite halo is not 
strongly affected by the tides from the host halo. However, 



owing to their rather complicated configuration, they could 
not discern t he effect of the reson ant torque, although some 
authors (e.g. lHavashi et al1l2003f ) has noticed that analytic 
formulae underestimate the mass loss. 

Although we have improved our understanding of the 
detailed physical processes responsible for satellite disrup- 
tion, some issues remain. By separating the heating mech- 
anisms into two distinct regimes, we achieved an improved 
understanding of resonant dynamics for an eccentric orbit. 
However, we have not compared this approximation with 
a comprehensive perturbation theory calculation. The res- 
onant heating of satellites on eccentric orbits is similar to 
heating by other subhaloes. This interaction is an impor- 
tant source of satellite evolution in addition to the inter- 
action with the smooth host halo. In addition, the initially 
spherical satellite is deformed during disruption owing to the 
host halo's tidal field and we have not yet accounted for this 
deformation. A comprehensive treatment of this deforma- 
tion might be necessary to understand satellite disruption 
in detail. Lastly, we also need a better understanding of the 
non-linear processes that occur during satellite evolution. 

In this study, we characterised the linear processes; un- 
derstanding the detailed consequences of the non-linear pro- 
cesses is a daunting future task. We should then finally be 
able to fully constrain the satellite disruption mechanism, 
which is an essential ingredient of galaxy formation and evo- 
lution. 
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APPENDIX A: A SATELLITE'S EFFECTIVE 
POTENTIAL IN A GENERAL 

GRAVITATIONAL FIELD AND A UNIFORMLY 
ROTATING FRAME 

The effective potential of a spherical satellite as function of 
a satellite's radius (r) is determined by its self-gravity, the 
external potential, and the centrifugal force. We set r — 
R — Ro where R is the location relative to the host halo's 
centre and Ro is the location of the satellite's centre relative 
to the host halo's centre. The acceleration in the satellite's 
frame, r, is difference between the acceleration, R, and the 
effective acceleration of the satellite, Ro, in the host frame. 
The quantities R and Ro are 



R = -V$tot(R) -tlx {Six R), 
Ro = -V$ to *(-Ro) - ft x {ft x Ro), 



(Al) 
(A2) 



where $ tot (Ji) is the total potential on the particle and ft 
is angular velocity of the satellite. The quantity & t ot{R) = 
® host {R) + $sat{R), where & ho3t {R) is the host halo po- 
tential and & sa t{R) is the potential of the satellite. Using 
r = R — Ro, the equation of motion in the satellite frame 
becomes 

r = -V$ sat (r) - V$ ho3 t{R) 

-V$ host (i? ) -nx(nxr). (A3) 

Because we assume a circular orbit, & sa t{Ro) is a 
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constant; V$ sa t(-R) -> V$ sa i(r) and V$t ot (Ro) - 
V$h ost (i?o). The second and third terms in equation ()A3 
can be further simplified using a Taylor expansion, 



V$ host (i?) - V$ host (i? ) 



d 2 $ host (R ) 
dR 2 



r + O [{r/Ro) 



(A4) 



The last term in equation ()A3[) is the centrifugal term. As- 
suming that the orbital plane is equatorial with Cl = z and 
r = (x,y,x), we have 

- S7 x (O x r) = -[(fi ■ r)n - (n • fl)r] = Q 2 {x,y,0).{A5) 

For generating initial conditions in §2]and our idealised mod- 
els in [J6] we will approximate equation (|A5|) with a spherical 
average 



- fl x (ft x r) » aff r 

where a is between and 1. 

Using this expansion, equation ()A3[1 becomes 



(A6) 



-V$ sat (r) 



d 2 $ hosf (R ] 
dR 2 



r + qS7 r 



= -V$««t(r) - [47rGp ftost (i?o) - (2 + a)fl 2 ]r. (A7) 

Finally, a satellite's effective potential as function of r using 
equation (|A7j) can be written as: 



$ e// (r) = $„t(r) 

+ [47rGp host ( 7?o ) - (2 -!-,>)<>-] < ( ) . (A.sj 



REFERENCES 

Besla G., Kallivayalil N., Hernquist L., Robertson B., Cox 
T. J., van der Marel R. P., Alcock C, 2007, ApJ, 668, 949 

Binney J., Tremaine S., 1987, Galactic dynamics. Prince- 
ton, NJ, Princeton University Press, 1987 

Boylan-Kolchin M., Ma C.-R, 2007, MNRAS, 374, 1227 

Bullock J. S., Kravtsov A. V., Weinberg D. H., 2001, ApJ, 
548, 33 

Chernoff D. F., Weinberg M. D., 1990, ApJ, 351, 121 
Choi J.-H., Weinberg M. D., Katz N., 2007, MNRAS, 381, 
987 

Clutton-Brock M., 1972, Ap&SS, 16, 101 
Clutton-Brock M., 1973, Ap&SS, 23, 55 
Connors T. W., Kawata D., Gibson B. K., 2006, MNRAS, 
371, 108 

De Lucia G., Kauffmann G., Springel V., White S. D. M., 
Lanzoni B., Stoehr F., Tormen G., Yoshida N., 2004, MN- 
RAS, 348, 333 
Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859 
Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., 

Potter D., Stadel J., 2008, Nature, 454, 735 
Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535 
Gao L., White S. D. M., Jenkins A., Stoehr F., Springel 

V., 2004, MNRAS, 355, 819 
Gardiner L. T., Noguchi M., 1996, MNRAS, 278, 191 
Gardiner L. T., Sawa T., Fujimoto M., 1994, MNRAS, 266, 
567 

Ghigna S., Moore B., Governato F., Lake G., Quinn T., 
Stadel J., 1998, MNRAS, 300, 146 



Ghigna S., Moore B., Governato F., Lake G., Quinn T., 

Stadel J., 2000, ApJ, 544, 616 
Gnedin O. Y., Hernquist L., Ostriker J. P., 1999, ApJ, 514, 

109 

Gnedin O. Y., Ostriker J. P., 1999, ApJ, 513, 626 
Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn 

T., 2003, ApJ, 584, 541 
Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375 
Kallivayalil N., van der Marel R. P., Alcock C, Axelrod T., 

Cook K. H., Drake A. J., Geha M., 2006, ApJ, 638, 772 
Kazantzidis S., Mayer L., Mastropietro C, Diemand J., 

Stadel J., Moore B., 2004, ApJ, 608, 663 
Kim S., Staveley-Smith L., Dopita M. A., Freeman K. C, 

Sault R. J., Kesteven M. J., McConnell D., 1998, ApJ, 

503, 674 

Klypin A., Gottlober S., Kravtsov A. V., Khokhlov A. M., 

1999, ApJ, 516, 530 
Klypin A., Zhao H., Somerville R. S., 2002, ApJ, 573, 597 
Kunkel W. E., Demers S., Irwin M. J., Albert L., 1997, 

ApJL, 488, L129+ 
Majewski S. R. et al., AJ, 128, 245 

Mastropietro C, Moore B., Mayer L., Wadsley J., Stadel 

J., 2005, MNRAS, 363, 509 
Mathewson D. S., Schwarz M. P., Murray J. D., 1977, 

ApJL, 217, L5 
Moore B., Davis M., 1994, MNRAS, 270, 209 
Moore B., Katz N., Lake G., 1996, ApJ, 457, 455 
Murai T., Fujimoto M., 1980, PASJ, 32, 581 
Murali C, Katz N., Hernquist L., Weinberg D. H., Dave 

R., 2002, ApJ, 571, 1 
Murray C. D., Dermott S. F., 1999, Solar system dynamics. 

Solar system dynamics by Murray, C. D., 1999 
Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 

493 

Oguri M., Lee J., 2004, MNRAS, 355, 120 

Penarrubia J., Navarro J. F., McConnachie A. W., 2008, 

ApJ, 673, 226 
Putman M. E. et al. 1998, Nature, 394, 752 
Quilis V., Moore B., Bower R., 2000, Science, 288, 1617 
Read J. I., Wilkinson M. L, Evans N. W., Gilmore G., 

Kleyna J. T., 2006, MNRAS, 367, 387 
Spitzer L., 1987, Dynamical evolution of globular clusters. 

Princeton, NJ, Princeton University Press, 1987 
Springel V., Wang J., Vogelsberger M., Ludlow A., Jenkins 

A., Helmi A., Navarro J. F., Frenk C. S., White S. D. M., 

2008, MNRAS, 391, 1685 
Stoehr F., 2006, MNRAS, 365, 147 

Stoehr F., White S. D. M., Tormen G., Springel V., 2002, 

MNRAS, 335, L84 
Taylor J. E., Babul A., 2001, ApJ, 559, 716 
Taylor J. E., Babul A., 2004, MNRAS, 348, 811 
Taylor J. E., Babul A., 2005a, MNRAS, 364, 515 
Taylor J. E., Babul A., 2005b, MNRAS, 364, 535 
Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729 
van der Marel R. P., Alves D. R., Hardy E., Suntzeff N. B., 

2002, AJ, 124, 2639 
Weinberg M. D., 1986, ApJ, 300, 93 
Weinberg M. D., 1989, MNRAS, 239, 549 
Weinberg M. D., 1994a, AJ, 108, 1398 
Weinberg M. D., 1994b, AJ, 108, 1403 
Weinberg M. D., 1999 AJ, 117, 629 
Weinberg M. D., 2004, |astro-ph/0404169l 



© 0000 RAS, MNRAS 000, 000-000 



20 Choi, Weinberg, & Katz 

Weinberg M. D., Katz N., 2007a, MNRAS, 375, 425 
Weinberg M. D., Katz N., 2007b, MNRAS, 375, 460 
Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., 
Wechsler R. H., 2005, ApJ, 624, 505 



© 0000 RAS, MNRAS 000, 000-000 



0.0016 



0.0014 



0.0012 



0.001 



<B 0.0008 



0.0006 



0.0004 



0.0002 




